rm(list = ls())

library(R2jags)

###
###  Data
###

load("data.RData")

data$Q5.1 <- factor(data$Q5.1, levels = c(1:5, 98, 99)) 
data$Q5.1A <- factor(data$Q5.1A, levels = c(1:5, 98, 99)) 
data$Q5.1B <- factor(data$Q5.1B, levels = c(1:5, 98, 99)) 
data$Q5.2 <- factor(data$Q5.2, levels = c(1:5, 98, 99)) 
data$Q5.2A <- factor(data$Q5.2A, levels = c(1:5, 98, 99)) 
data$Q5.2B <- factor(data$Q5.2B, levels = c(1:5, 98, 99)) 
data$Q5.3 <- factor(data$Q5.3, levels = c(1:5, 98, 99)) 
data$Q5.3A <- factor(data$Q5.3A, levels = c(1:5, 98, 99)) 
data$Q5.3B <- factor(data$Q5.3B, levels = c(1:5, 98, 99)) 
data$Q5.4 <- factor(data$Q5.4, levels = c(1:5, 98, 99)) 
data$Q5.4A <- factor(data$Q5.4A, levels = c(1:5, 98, 99)) 
data$Q5.4B <- factor(data$Q5.4B, levels = c(1:5, 98, 99)) 

###
###  Descriptive Plots by Provinces
###

yaxis <- c("Control", "ISAF", "Taliban")
col.seq <- c(grey(seq(0.3^2.2, 0.9^2.2, length = 5)^(1/2.2)), rep(rgb(0,0,0,0), 2)) 

n.samp <- c(nrow(data), rep(NA, 5))

pdf("figs/Fig1-EndorseDescriptiveProvinces.pdf", height = 5.5, width = 6.6)
par(mfrow = c(6,4), las = 1, mar=c(0,0.5,0.5,0.25), oma=c(4,6,2,0.5), cex = 0.55)

barplot(t(rbind(matrix(prop.table(table(data$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = yaxis, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = yaxis, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

par(xpd = NA)
lines(c(-3.8, 1.2), rep(-0.125, 2), lty = "dashed")

for (i in 1:length(levels(data$M5))) {

  dat <- data[data$M5 == levels(data$M5)[i], ]

  n.samp[i+1] <- nrow(dat)

  barplot(t(rbind(matrix(prop.table(table(dat$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = yaxis, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = yaxis, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
}

mtext("Direct Elections", outer = TRUE, at = 0.125, cex = 0.75, line = -0.25)
mtext("Prison Reform", outer = TRUE, at = 0.375, cex = 0.75, line = -0.25)
mtext("Independent Election\n Commission", outer = TRUE, at = 0.625, cex = 0.75, line = -0.5)
mtext("Anti-Corruption\n Reform", outer = TRUE, at = 0.88, cex = 0.75, line = -0.5)

par(las = 0)
labels <- c("Overall", levels(data$M5))
for (i in 1:length(labels)) {
  mtext(labels[i], outer = TRUE, side = 2, line = 4.5, at = 0.9 - 0.165*(i-1), cex = 0.75)
  mtext(paste("(N = ", n.samp[i],")", sep=""), outer = TRUE, side = 2, line = 3.5, at = 0.9 - 0.165*(i-1), cex = 0.5)
}

par(xpd=NA)
legend(-1, -0.25, c(" Strongly\n agree", " Agree", " Indifferent", " Disgree",
                       " Strongly\n disagree", " Don't Know", " Refused"),
       density = c(rep(-1, 5), rep(40, 2)), angle = c(rep(-45, 6), 45),
       fill = c(col.seq[1:5], rep(rgb(.31, .31, .31), 2)), horiz=T, xjust=0.5,
       x.intersp = 0.25, bty = "n")

dev.off()



###
###  Descriptive Plots for Helmand
###

i <- 1

pdf(paste("figs/Fig8-EndorseDescriptive", levels(data$M5)[i], ".pdf", sep = ""),
    height = 4.75, width = 6.6)

dat <- data[data$M5 == levels(data$M5)[i], ]
n.samp <- rep(NA, length(unique(dat$DIS)))

par(mfrow = c(5,4), las = 1, mar=c(0,0.5,0.5,0.25), oma=c(4,6,2,0.5), cex = 0.55)
print(length(unique(dat$DIS)))

for (j in 1:length(unique(dat$DIS))) {
  dat.sub <- dat[dat$DIS == unique(dat$DIS)[j], ]
  n.samp[j] <- nrow(dat.sub)
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = yaxis, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.1)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.1A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.1B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = yaxis, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.2)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.2A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.2B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.3)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.3A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.3B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat.sub$Q5.4)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.4A)[c(5:1, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat.sub$Q5.4B)[c(5:1, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
 
  
}


mtext("Direct Elections", outer = TRUE, at = 0.125, cex = 0.75, line = -0.25)
mtext("Prison Reform", outer = TRUE, at = 0.375, cex = 0.75, line = -0.25)
mtext("Independent Election\n Commission", outer = TRUE, at = 0.625, cex = 0.75, line = -0.5)
mtext("Anti-Corruption\n Reform", outer = TRUE, at = 0.88, cex = 0.75, line = -0.5)

par(las = 0)
for (j in 1:length(unique(dat$DIS))) {
  mtext(paste(unique(dat$DIS)[j]), outer = TRUE, side = 2, line = 4.5, at = 0.89 - 0.2*(j-1), cex = 0.75)
  mtext(paste("(N = ", n.samp[j],")", sep=""), outer = TRUE, side = 2, line = 3.5, at = 0.89 - 0.2*(j-1), cex = 0.5)
}

mtext(paste(levels(data$M5)[i], "Province"), outer = TRUE, side = 2, line = 6, at = 0.49, cex = 1)

par(xpd=NA)
legend(-1, -0.25, c(" Strongly\n agree", " Agree", " Indifferent", " Disgree",
                    " Strongly\n disagree", " Don't Know", " Refused"),
       density = c(rep(-1, 5), rep(40, 2)), angle = c(rep(-45, 6), 45),
       fill = c(col.seq[1:5], rep(rgb(.31, .31, .31), 2)), horiz=T, xjust=0.5,
       x.intersp = 0.25, bty = "n")

dev.off()

## #########
## Figure 2

pdf(file = "figs/Fig2-DemographicsHistograms.pdf", height = 2, width = 7.5)

par(mfcol=c(1,4), las = 1, mar=c(1.5,1.5,2,2), oma=c(2,1.5,2,0), cex = 0.525, cex.axis = 1)

age.cat <- data$Q1.1 <= 18
age.cat[data$Q1.1 > 18 & data$Q1.1 <= 25] <- 2
age.cat[data$Q1.1 > 25 & data$Q1.1 <= 35] <- 3
age.cat[data$Q1.1 > 35 & data$Q1.1 <= 45] <- 4
age.cat[data$Q1.1 > 45 & data$Q1.1 <= 60] <- 5
age.cat[data$Q1.1 > 60] <- 6

barplot(table(age.cat), names.arg = c("<18", "25", "35","45","60", ">60"), ylim = c(0, 1400))
abline(v = 3.1, lty = "dashed")
text(4.2, 1200, "Median")

educ <- data$Q1.2
educ[data$Q1.2 >= 12] <- 12

barplot(table(educ), axes = F, names.arg = NA, ylim = c(0, 1400))
abline(v = 1*1.2+.7, lty = "dashed")

axis(1, at = c(0, 2, 4, 6, 8, 10, 12)*1.2+0.7, labels = c("0", "2", "4", "6", "8", "10", "12+"),
     xpd = TRUE, lwd = 0)
axis(2)
     
madrassa <- data$Q1.4
madrassa[is.na(data$Q1.4)] <- 0

data$Q1.3[data$Q1.3 > 1] <- NA
data$Q1.4[data$Q1.3 == 0] <- 0
madr <- data$Q1.4
madr[data$Q1.4 > 6] <- 7

barplot(table(madr), axes = F, names.arg = NA, ylim = c(0, 1400))
axis(1, at = (1:8)*1.2-.5, labels = c(0:6, "7+"), lwd = 0)
axis(2)
abline(v = 2*1.2 - .5, lty = "dashed")

barplot(table(data$income), ylim = c(0, 1400), names.arg = NA, axes = F)
axis(1, at = (1:5)*1.2-.5, labels = c("<\n2k", "2-\n10k", "10-\n20k", "20-\n30k", ">\n30k"),
     mgp = c(3, 2, 0), lwd = 0)
axis(2)
abline(v = 2*1.2-.5, lty = "dashed")

mtext("Age", side = 3, outer = TRUE, at = 0.125, line = -1, cex = 0.7)
mtext("Years of Education", side = 3, outer = TRUE, at = 0.375, line = -1, cex = 0.7)
mtext("Years of Madrassa\n Schooling", side = 3, outer = TRUE, at = 0.625, line = -1, cex = 0.7)
mtext("Income (Afghanis)", side = 3, outer = TRUE, at = 0.875, line = -1, cex = 0.7)

dev.off()


## ####
## Figure 11 Balance

pdf(file = "figs/Fig9-BalanceSample.pdf", height = 2, width = 7.125)
par(mfcol = c(1,4), las = 1, mar=c(1.5,1.5,2,2), oma=c(2,3,2,0), cex = 0.525)
plot(density(villages.balance$altitude[villages.balance$village.surveyed==1]/1000, na.rm = T), xlab = "Altitude (hundreds meters)", ylab = "Density", main = "", ylim = c(0,1), axes = F, frame = TRUE)
lines(density(villages.balance$altitude[villages.balance$village.surveyed==0 & villages.balance$district.surveyed == 1]/1000, na.rm = T), col = "darkgray", lwd = 1)
axis(2)
axis(1, at = c(0,1,2,3))


plot(density(villages.balance$population[villages.balance$village.surveyed==1]/100, na.rm = T), xlab = "Population count (hundreds)", ylab = "Density",main = "", ylim = c(0, .15))
lines(density(villages.balance$population[villages.balance$village.surveyed==0 & villages.balance$district.surveyed == 1]/100, na.rm = T), col = "darkgray", lwd = 1)

plot(density(villages.balance$violence.count.ISAF.5km[villages.balance$village.surveyed==1], na.rm = T), xlab = "Event count (within 5km)", ylab = "Density", ylim=c(0,.15), main = "")
lines(density(villages.balance$violence.count.ISAF.5km[villages.balance$village.surveyed==0 & villages.balance$district.surveyed == 1], na.rm = T), col = "darkgray", lwd = 1)

plot(density(villages.balance$violence.count.taliban.5km[villages.balance$village.surveyed==1], na.rm = T), xlab = "Event count (within 5km)", ylab = "Density", ylim=c(0,0.025), main = "")
lines(density(villages.balance$violence.count.taliban.5km[villages.balance$village.surveyed==0 & villages.balance$district.surveyed == 1], na.rm = T), col = "darkgray", lwd = 1)

legend("topright", c("Surveyed Village","Non Surveyed in\nSampled Districts"),
       lty = c("solid","solid"), col = c("black", "darkgray"),
       bty = "n")

mtext("Altitude (km)", side = 3, outer = TRUE, at = 0.125, line = -1, cex = 0.7)
mtext("Population Size (hundreds)", side = 3, outer = TRUE, at = 0.375, line = -1, cex = 0.7)
mtext("Violent Events\nInitiated by ISAF", side = 3, outer = TRUE, at = 0.625, line = -1.5, cex = 0.7)
mtext("Violent Events\nInitiated by the Taliban", side = 3, outer = TRUE, at = 0.875, line = -1.5, cex = 0.7)

dev.off()


## ######
## Fig 12

data$collapsed.ISAF.indiv <- factor(data$violent.exp.ISAF * (data$encounter.ISAF + 1), levels = c(2,1,0), labels = c( "Approach", "No approach", "No harm"))
data$collapsed.taliban.indiv <- factor(data$violent.exp.taliban * (data$encounter.taliban + 1), levels = c(2,1,0), labels = c( "Approach", "No approach", "No harm"))

data$collapsed.ISAF.area <- factor(data$violent.exp.area.ISAF * (data$encounter.area.ISAF + 1), levels = c(2,1,0), labels = c( "Harm and\nApproach", "Harm and\nno approach", "No harm"))
data$collapsed.taliban.area <- factor(data$violent.exp.area.taliban * (data$encounter.area.taliban + 1), levels = c(2,1,0), labels = c( "Harm and\nApproach", "Harm and\nno approach", "No harm"))

data$collapsed.ISAF.area <- factor(data$violent.exp.area.ISAF * (data$encounter.area.ISAF + 1), levels = c(2,1,0), labels = c( "Harm and\nApproach", "Harm and\nno approach", "No harm"))
data$collapsed.taliban.area <- factor(data$violent.exp.area.taliban * (data$encounter.area.taliban + 1), levels = c(2,1,0), labels = c( "Harm and\nApproach", "Harm and\nno approach", "No harm"))

library(vcd)
vnames <- list(set_varnames = c(violent.exp.ISAF.fac="Harmed by ISAF",
                 violent.exp.taliban.fac="Harmed by Taliban",
                 encounter.taliban.fac="Approached by Taliban",
                 encounter.ISAF.fac = "Approached by ISAF"))
vnames <- list(set_varnames = c(violent.exp.ISAF.fac="",
                 violent.exp.taliban.fac="",
                 encounter.taliban.fac="",
                 encounter.ISAF.fac = ""))

show.prop <- .07

value.labels <- c( "Harmed and\nApproached", "Harmed and\nNot Approached", "Not Harmed")

pdf(file = "figs/Fig10A-ViolenceMosaicIndiv.pdf", width = 12.5, height = 15.3)

par(cex = .8)

vnames <- list(set_varnames = c(collapsed.ISAF.indiv="ISAF",
                collapsed.taliban.indiv ="Taliban"))

pushViewport(viewport(layout = grid.layout(nrow = 3, ncol = 2, just = c("center", "top"), widths = c(unit(6, "inches"), unit(4, "inches")), heights = c(unit(6.3, "inches"), unit(4.75, "inches"), unit(4.75, "inches")))))

pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1))

ISAFIndiv <- factor(data$violent.exp.ISAF * (data$encounter.ISAF + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoapproach", "NoHarm"))
TalibanIndiv <- factor(data$violent.exp.taliban * (data$encounter.taliban + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoApproach", "NoHarm"))

tab <- table(ISAFIndiv, TalibanIndiv)

mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanIndiv", highlighting_fill = grey.colors,
                       labeling_args = list(rot_labels = c(left = 0, top = 25), offset_labels = c(left = 0.3, top = 1), offset_varnames = c(left = 5.5, top = 2.2),
                       set_varnames = c(ISAFIndiv = "ISAF", TalibanIndiv ="Taliban"), 
                       set_labels = list(ISAFIndiv = value.labels,
                         TalibanIndiv = value.labels),
                       gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                         gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                         gp_main = gpar(fontface = "bold"),
                         just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, ##keep_aspect_ratio = TRUE,
                     main = paste("Overall (N = ", nrow(data), ")", sep = ""),
                       margins = c(left = 10, top = 2, bottom = -4.5, right = 2), pop = FALSE)

tabDisplay <- round(ifelse( tab / apply(tab, 1, sum) > .1 & tab / apply(tab, 2, sum) > .1, tab, NA) / sum(tab), 2)
labeling_cells(text = tabDisplay, margin = 0)(tab)

popViewport()

popViewport()

for(prov.id in 1:5) {
  if(prov.id == 1)
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 2))
  else if (prov.id == 2)
    pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 1))
  else if (prov.id == 3)
    pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
  else if (prov.id == 4)
    pushViewport(viewport(layout.pos.row = 3, layout.pos.col = 1))
  else if (prov.id == 5)
    pushViewport(viewport(layout.pos.row = 3, layout.pos.col = 2))

  print(prov.id)
  
  data.subset <- data[data$province.id == prov.id,]

  ISAFIndiv <- factor(data.subset$violent.exp.ISAF * (data.subset$encounter.ISAF + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoapproach", "NoHarm"))
  TalibanIndiv <- factor(data.subset$violent.exp.taliban * (data.subset$encounter.taliban + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoApproach", "NoHarm"))

  tab <- table(ISAFIndiv, TalibanIndiv)

  if (prov.id == 1)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanIndiv", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFIndiv = "ISAF", TalibanIndiv ="Taliban"), 
                           set_labels = list(ISAFIndiv = value.labels,
                             TalibanIndiv = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                         gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = 2.7, top = 3.3, bottom = -0.6, right = 2.55), pop = FALSE)
  
  else if (prov.id == 2 | prov.id == 4)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanIndiv", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFIndiv = "ISAF", TalibanIndiv ="Taliban"), 
                           set_labels = list(ISAFIndiv = value.labels,
                             TalibanIndiv = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                           gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = -2, top = 1, bottom = 1, right = -10), pop = FALSE)
  else if (prov.id == 3 | prov.id == 5)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanIndiv", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFIndiv = "ISAF", TalibanIndiv ="Taliban"), 
                           set_labels = list(ISAFIndiv = value.labels,
                             TalibanIndiv = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                           gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = 0, top = 1, bottom = 1, right = -0.5), pop = FALSE)

  
tabDisplay <- round(ifelse( tab / apply(tab, 1, sum) > .1 & tab / apply(tab, 2, sum) > .1, tab, NA) / sum(tab), 2)
labeling_cells(text = tabDisplay, margin = 0)(tab)

  popViewport()
  popViewport()
  
}

dev.off()
##

pdf(file = "figs/Fig10B-ViolenceMosaicArea.pdf", width = 12.5, height = 15.9)

pushViewport(viewport(layout = grid.layout(nrow = 3, ncol = 2, just = c("center", "top"), widths = c(unit(6, "inches"), unit(4, "inches")), heights = c(unit(6.5, "inches"), unit(4.95, "inches"), unit(4.95, "inches")))))##, heights = c(1,0.5,0.5))))

pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1))

ISAFArea <- factor(data$violent.exp.area.ISAF * (data$encounter.area.ISAF + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoapproach", "NoHarm"))
TalibanArea <- factor(data$violent.exp.area.taliban * (data$encounter.area.taliban + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoApproach", "NoHarm"))

tab <- table(ISAFArea, TalibanArea)

mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanArea", highlighting_fill = grey.colors,
                       labeling_args = list(rot_labels = c(left = 0, top = 25), offset_labels = c(left = 0.3, top = 1), offset_varnames = c(left = 5.5, top = 2.2),
                       set_varnames = c(ISAFArea = "ISAF", TalibanArea ="Taliban"), 
                       set_labels = list(ISAFArea = value.labels,
                         TalibanArea = value.labels),
                       gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                         gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                         gp_main = gpar(fontface = "bold"),
                         just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, ##keep_aspect_ratio = TRUE,
                     main = paste("Overall (N = ", nrow(data), ")", sep = ""),
                       margins = c(left = 10, top = 4, bottom = -0.3, right = 2), pop = FALSE)

tabDisplay <- round(ifelse( tab / apply(tab, 1, sum) > .1 & tab / apply(tab, 2, sum) > .1, tab, NA) / sum(tab), 2)
labeling_cells(text = tabDisplay, margin = 0)(tab)

popViewport()

popViewport()

for(prov.id in 1:5) {
  
  if(prov.id == 1)
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 2))
  else if (prov.id == 2)
    pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 1))
  else if (prov.id == 3)
    pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
  else if (prov.id == 4)
    pushViewport(viewport(layout.pos.row = 3, layout.pos.col = 1))
  else if (prov.id == 5)
    pushViewport(viewport(layout.pos.row = 3, layout.pos.col = 2))
  
  print(prov.id)
  
  data.subset <- data[data$province.id == prov.id,]

  ISAFArea <- factor(data.subset$violent.exp.area.ISAF * (data.subset$encounter.area.ISAF + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoapproach", "NoHarm"))
  TalibanArea <- factor(data.subset$violent.exp.area.taliban * (data.subset$encounter.area.taliban + 1), levels = c(2,1,0), labels = c( "HarmApproach", "HarmNoApproach", "NoHarm"))

  tab <- table(ISAFArea, TalibanArea)

   
  if (prov.id == 1)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanArea", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFArea = "ISAF", TalibanArea ="Taliban"), 
                           set_labels = list(ISAFArea = value.labels,
                             TalibanArea = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                         gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = 2.7, top = 3.3, bottom = -0.6, right = 2.55), pop = FALSE)
  
  else if (prov.id == 2 | prov.id == 4)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanArea", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFArea = "ISAF", TalibanArea ="Taliban"), 
                           set_labels = list(ISAFArea = value.labels,
                             TalibanArea = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                           gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = -2, top = 1, bottom = 1, right = -10), pop = FALSE)
  else if (prov.id == 3 | prov.id == 5)
    mosaic.obj <- mosaic(tab, spacing = spacing_equal(sp = unit(0.15, "inches")), highlighting = "TalibanArea", highlighting_fill = grey.colors,
                         labeling_args = list(rot_labels = c(left = 0, top = 35), offset_labels = c(left = 500, top = 500), offset_varnames = c(left = 500, top = 500),
                           set_varnames = c(ISAFArea = "ISAF", TalibanArea ="Taliban"), 
                           set_labels = list(ISAFArea = value.labels,
                             TalibanArea = value.labels),
                           gp_labels = gpar(lineheight = 0.8, fontsize = 16),
                           gp_varnames = gpar(fontsize = 17, fontface = "bold"),
                           gp_main = gpar(fontface = "bold"),
                           just_labels = c(left = "right", top = "left"), value_type = "observed"), newpage = FALSE, cex = .8, keep_aspect_ratio = TRUE,
                         main = paste(unique(data$M5[data$province.id == prov.id]), " (N = ", sum(data$province.id == prov.id), ")", sep=""),
                         margins = c(left = 0, top = 1, bottom = 1, right = -0.5), pop = FALSE)

  tabDisplay <- round(ifelse( tab / apply(tab, 1, sum) > .1 & tab / apply(tab, 2, sum) > .1, tab, NA) / sum(tab), 2)
  
  labeling_cells(text = tabDisplay, margin = 0)(tab)
  
  popViewport()
  popViewport()
  
}

dev.off()

## #####
## create demographic - support plot

load("results/endorseResultsNoCovariates.RData")

attach.jags(fit)

s.mean <- apply(s, c(1,2), mean)
sd.x.matrix <-  matrix(sd.x, byrow = F, nrow = nrow(s.mean), ncol = ncol(s.mean))
s.mean.std <- s.mean / sd.x.matrix

s.indiv <- apply(s.mean.std, 2, mean)
s.indiv.sd <- apply(s.mean.std, 2, sd)

data$s.indiv <- s.indiv
data$s.indiv.weights <- 1/s.indiv.sd^2
data$endorser.ISAF <- (data$endorser == 2)*1
data$endorser.taliban <- (data$endorser == 3)*1

data$violent.exp.int <- data$violent.exp.taliban * data$violent.exp.ISAF
data.isaf <- data[data$endorser.ISAF==1,]
data.taliban <- data[data$endorser.taliban==1,]

pdf("figs/Fig4-Demographics.pdf", height = 4, width = 7.125)
par(mfcol = c(2,4), las = 1, mar=c(1.5,1.5,2,2), oma=c(2,3,2,0), cex = 0.525)

## AGE: older => less support for Taliban, no clear relationship for ISAF
plot(data.taliban$s.indiv ~ data.taliban$Q1.1, pch = ".", axes = FALSE, frame = TRUE,
     cex = 2, xlim = c(10, 80), ylim = c(-3,2))
lines(lowess(data.taliban$s.indiv ~ data.taliban$Q1.1), lwd = 2)
abline(0, 0, lty = "dashed")
axis(1, at = seq(from = 10, to = 80, by = 10))     
axis(2)

plot(data.isaf$s.indiv ~ data.isaf$Q1.1, pch = ".", cex = 2, axes = FALSE, frame = TRUE,
     xlim = c(10, 80), ylim = c(-3,2))
lines(lowess(data.isaf$s.indiv ~ data.isaf$Q1.1), lwd = 2)
abline(0, 0, lty = "dashed")
axis(1, at = seq(from = 10, to = 80, by = 10))     
axis(2)

## EDUCATION: More educated => less support for Taliban, more support for ISAF
## very clear, almost linear pattern
educ.taliban <- data.taliban$Q1.2
educ.taliban[data.taliban$Q1.2 >= 12] <- 12
educ.isaf <- data.isaf$Q1.2
educ.isaf[data.isaf$Q1.2 >= 12] <- 12

boxplot(data.taliban$s.indiv ~ educ.taliban, varwidth = TRUE, ylim = c(-3,2),
        axes = FALSE, frame = TRUE)
abline(0, 0, lty = "dashed")
axis(1, at = c(0, 2, 4, 6, 8, 10, 12)+1, labels = c("0", "2", "4", "6", "8", "10", "12+"),
     xpd = TRUE)
axis(2)

boxplot(data.isaf$s.indiv ~ educ.isaf, varwidth = TRUE, ylim = c(-3,2), axes = FALSE,
        frame = TRUE)
abline(0, 0, lty = "dashed")
axis(1, at = c(0, 2, 4, 6, 8, 10, 12)+1, labels = c("0", "2", "4", "6", "8", "10", "12+"),
     xpd = TRUE)
axis(2)

## MADRASSA SCHOOLING: Yes => more support for Taliban, less support for ISAF
## More schooling => more support for Taliban, no clear relationship for ISAF
data.isaf$Q1.3[data.isaf$Q1.3 > 1] <- NA
data.taliban$Q1.3[data.taliban$Q1.3 > 1] <- NA
data.isaf$Q1.4[data.isaf$Q1.3 == 0] <- 0
data.taliban$Q1.4[data.taliban$Q1.3 == 0] <- 0
madr.taliban <- data.taliban$Q1.4
madr.taliban[data.taliban$Q1.4 > 6] <- 7
madr.isaf <- data.isaf$Q1.4
madr.isaf[data.isaf$Q1.4 > 6] <- 7

boxplot(data.taliban$s.indiv ~ madr.taliban, varwidth = TRUE, ylim = c(-3,2), axes = FALSE,
        frame = TRUE)
abline(0, 0, lty = "dashed")
axis(1, at = 1:8, labels = c(0:6, "7+"))
axis(2)

boxplot(data.isaf$s.indiv ~ madr.isaf, varwidth = TRUE, ylim = c(-3,2), axes = FALSE,
        frame = TRUE)
abline(0, 0, lty = "dashed")
axis(1, at = 1:8, labels = c(0:6, "7+"))
axis(2)

## INCOME: Richer => less support for Taliban, more support for ISAF
## again very clear, almost linear pattern (except perhaps the last category)
##data.isaf$Q1.12[data.isaf$Q1.12 > 5] <- NA
##data.taliban$Q1.12[data.taliban$Q1.12 > 5] <- NA
boxplot(data.taliban$s.indiv ~ data.taliban$income, varwidth = TRUE, ylim = c(-3,2),
        frame = TRUE, axes = FALSE)
abline(0, 0, lty = "dashed")
axis(1, at = 1:5, labels = c("<\n2k", "2-\n10k", "10-\n20k", "20-\n30k", ">\n30k"),
     mgp = c(3, 2, 0))
axis(2)

boxplot(data.isaf$s.indiv ~ data.isaf$income, varwidth = TRUE, ylim = c(-3,2), frame = TRUE,
        axes = FALSE)
abline(0, 0, lty = "dashed")
axis(1, at = 1:5, labels = c("<\n2k", "2-\n10k", "10-\n20k", "20-\n30k", ">\n30k"),
     mgp = c(3, 2, 0))
axis(2)

mtext("Age", side = 3, outer = TRUE, at = 0.125, line = -1, cex = 0.7)
mtext("Years of Education", side = 3, outer = TRUE, at = 0.375, line = -1, cex = 0.7)
mtext("Years of Madrassa\n Schooling", side = 3, outer = TRUE, at = 0.625, line = -1.5, cex = 0.7)
mtext("Income (Afghanis)", side = 3, outer = TRUE, at = 0.875, line = -1, cex = 0.7)

mtext("Support for Taliban", side = 2, outer = TRUE, las = 3, at = 0.75, line = 1.75, cex = 0.7)
mtext("Support for ISAF", side = 2, outer = TRUE, las = 3, at = 0.25, line = 1.75, cex = 0.7)
dev.off()


## ######
## Set up regression findings and calculate quantities of interest

rm(list=ls())

library(R2jags)
library(foreign)

for(model.iter in 1:3) {
  
if(model.iter != 2)
  violence.exp <- "individual"
  
if(model.iter == 2)
  violence.exp <- "area"
  
if(model.iter == 3)
  pro.taliban.model <- TRUE

if(exists("pro.taliban.model"))
  load("endorseSetupProTaliban.RData")
else
  load("endorseSetup.RData")

if (exists("no.covariate.model") == TRUE) {
  file.name <- "NoCovariates"
  violence.exp <- "individual"
} else if (exists("harm.only") == TRUE) {
  indivcov <- indivcov[, 1:4]
  file.name <- paste("HarmOnly.",violence.exp,".fit",sep="")
} else if (exists("prov.id") == TRUE) {
  violence.exp <- "individual"
  file.name <- paste("Prov.",prov.id,".fit",sep="")

  ## set up data for province model
  province.ind <- data$province.id == prov.id

  indivcov <- cbind(indivcov, indivcov[,1] * province.ind, indivcov[,3] * province.ind, indivcov[,5] * province.ind, indivcov[,7] * province.ind)
  colnames(indivcov)[17:20] <- paste("prov.",prov.id,c(".violent.taliban.int", ".violent.ISAF.int", ".encounter.taliban.int", ".encounter.ISAF.int"), sep = "")
  
} else if (exists("pro.taliban.model") == TRUE) {
  file.name <- paste("ProTal.",violence.exp,".fit",sep="")
  ## set up data for pro taliban model
  indivcov <- cbind(indivcov[,-16], indivcov[,1] * indivcov[,15], indivcov[,3] * indivcov[,15])
  colnames(indivcov)[16:17] <- c("pro.taliban.violent.taliban.int", "pro.taliban.violent.ISAF.int")
} else if (exists("pro.taliban.approach.model") == TRUE) {
  file.name <- paste("ProTalAp.",violence.exp,".fit",sep="")
  indivcov <- cbind(indivcov[,-16], indivcov[,1] * indivcov[,15], indivcov[,3] * indivcov[,15],
                            indivcov[,5] * indivcov[,15], indivcov[,7] * indivcov[,15])
  colnames(indivcov)[16:19] <- c("pro.taliban.violent.taliban.int", "pro.taliban.violent.ISAF.int",
                                 "pro.taliban.encounter.taliban.int", "pro.taliban.encounter.ISAF.int")
} else {
  if(exists("robustCheck") == TRUE) {
    file.name <- paste(robustCheck,".fit",sep="")
  } else {
    file.name <- paste(".",violence.exp,".fit",sep="")
  }
}

print(paste("file:",file.name))

if(exists("basic.model")) {
  indivcov <- indivcov[, 1:9]
  
  villcov <- as.matrix(rep(1, nrow(villcov)))
  
  distcov <- as.matrix(rep(1, nrow(distcov)))
}

for(iter in 1:3) {
  if(!exists("no.covariate.model")) {
    prefix <- "results/endorseResultsCovariates"
    if(!exists("pro.taliban.model"))
      prefix <- paste(prefix, sep = "")
    load(paste(prefix, file.name, ".", iter, ".RData", sep =""))
    if(iter == 1)
      fit.1 <- fit
    else if(iter == 2)
      fit.2 <- fit
    else if(iter == 3)
      fit.3 <- fit
  } else {
    load("results/endorseResultsNoCovariates.RData")
  }
}

attach.jags(fit.1)

theta.province.2.1 <- theta.province[,,2]
theta.province.3.1 <- theta.province[,,3]

theta.district.2.1 <- theta.district[,,2]
theta.district.3.1 <- theta.district[,,3]

lambda.village.2.1 <- lambda.village[,,2]
lambda.village.3.1 <- lambda.village[,,3]

lambda.cov1.1 <- lambda.cov1
theta.village.cov1.1 <- theta.village.cov1
theta.district.cov1.1 <- theta.district.cov1
alpha.1 <- alpha
beta.1 <- beta
delta.village.1 <- delta.village
delta.district.1 <- delta.district
delta.province.1 <- delta.province

sd.x.1 <- sd.x

detach.jags()

attach.jags(fit.2)

theta.province.2.2 <- theta.province[,,2]
theta.province.3.2 <- theta.province[,,3]

theta.district.2.2 <- theta.district[,,2]
theta.district.3.2 <- theta.district[,,3]

lambda.village.2.2 <- lambda.village[,,2]
lambda.village.3.2 <- lambda.village[,,3]

lambda.cov1.2 <- lambda.cov1
theta.village.cov1.2 <- theta.village.cov1
theta.district.cov1.2 <- theta.district.cov1
alpha.2 <- alpha
beta.2 <- beta
delta.village.2 <- delta.village
delta.district.2 <- delta.district
delta.province.2 <- delta.province

sd.x.2 <- sd.x

detach.jags()

attach.jags(fit.3)

theta.province.2.3 <- theta.province[,,2]
theta.province.3.3 <- theta.province[,,3]

theta.district.2.3 <- theta.district[,,2]
theta.district.3.3 <- theta.district[,,3]

lambda.village.2.3 <- lambda.village[,,2]
lambda.village.3.3 <- lambda.village[,,3]

lambda.cov1.3 <- lambda.cov1
theta.village.cov1.3 <- theta.village.cov1
theta.district.cov1.3 <- theta.district.cov1
alpha.3 <- alpha
beta.3 <- beta
delta.village.3 <- delta.village
delta.district.3 <- delta.district
delta.province.3 <- delta.province

sd.x.3 <- sd.x

detach.jags()

lambda.cov1.dim <- dim(lambda.cov1.1)
lambda.cov1.dim[1] <- dim(lambda.cov1.1)[1]/2 + dim(lambda.cov1.2)[1]/2 + dim(lambda.cov1.3)[1]/2 + 3
lambda.cov1 <- array(NA, dim = lambda.cov1.dim)
lambda.cov1[,1,] <- rbind(lambda.cov1.1[round(nrow(lambda.cov1.1)/2):nrow(lambda.cov1.1),1,],
                          lambda.cov1.2[round(nrow(lambda.cov1.2)/2):nrow(lambda.cov1.2),1,],
                          lambda.cov1.3[round(nrow(lambda.cov1.3)/2):nrow(lambda.cov1.3),1,])
lambda.cov1[,2,] <- rbind(lambda.cov1.1[round(nrow(lambda.cov1.1)/2):nrow(lambda.cov1.1),2,],
                          lambda.cov1.2[round(nrow(lambda.cov1.2)/2):nrow(lambda.cov1.2),2,],
                          lambda.cov1.3[round(nrow(lambda.cov1.3)/2):nrow(lambda.cov1.3),2,])
lambda.cov1[,3,] <- rbind(lambda.cov1.1[round(nrow(lambda.cov1.1)/2):nrow(lambda.cov1.1),3,],
                          lambda.cov1.2[round(nrow(lambda.cov1.2)/2):nrow(lambda.cov1.2),3,],
                          lambda.cov1.3[round(nrow(lambda.cov1.3)/2):nrow(lambda.cov1.3),3,])

theta.village.cov1.dim <- dim(theta.village.cov1.1)
theta.village.cov1.dim[1] <- dim(theta.village.cov1.1)[1]/2 + dim(theta.village.cov1.2)[1]/2 + dim(theta.village.cov1.3)[1]/2 + 3
theta.village.cov1 <- array(NA, dim = theta.village.cov1.dim)
theta.village.cov1[,1,] <- rbind(theta.village.cov1.1[round(nrow(theta.village.cov1.1)/2):nrow(theta.village.cov1.1),1,],
                          theta.village.cov1.2[round(nrow(theta.village.cov1.2)/2):nrow(theta.village.cov1.2),1,],
                          theta.village.cov1.3[round(nrow(theta.village.cov1.3)/2):nrow(theta.village.cov1.3),1,])
theta.village.cov1[,2,] <- rbind(theta.village.cov1.1[round(nrow(theta.village.cov1.1)/2):nrow(theta.village.cov1.1),2,],
                          theta.village.cov1.2[round(nrow(theta.village.cov1.2)/2):nrow(theta.village.cov1.2),2,],
                          theta.village.cov1.3[round(nrow(theta.village.cov1.3)/2):nrow(theta.village.cov1.3),2,])
theta.village.cov1[,3,] <- rbind(theta.village.cov1.1[round(nrow(theta.village.cov1.1)/2):nrow(theta.village.cov1.1),3,],
                          theta.village.cov1.2[round(nrow(theta.village.cov1.2)/2):nrow(theta.village.cov1.2),3,],
                          theta.village.cov1.3[round(nrow(theta.village.cov1.3)/2):nrow(theta.village.cov1.3),3,])

theta.district.cov1.dim <- dim(theta.district.cov1.1)
theta.district.cov1.dim[1] <- dim(theta.district.cov1.1)[1]/2 + dim(theta.district.cov1.2)[1]/2 + dim(theta.district.cov1.3)[1]/2 + 3
theta.district.cov1 <- array(NA, dim = theta.district.cov1.dim)
theta.district.cov1[,1,] <- rbind(theta.district.cov1.1[round(nrow(theta.district.cov1.1)/2):nrow(theta.district.cov1.1),1,],
                          theta.district.cov1.2[round(nrow(theta.district.cov1.2)/2):nrow(theta.district.cov1.2),1,],
                          theta.district.cov1.3[round(nrow(theta.district.cov1.3)/2):nrow(theta.district.cov1.3),1,])
theta.district.cov1[,2,] <- rbind(theta.district.cov1.1[round(nrow(theta.district.cov1.1)/2):nrow(theta.district.cov1.1),2,],
                          theta.district.cov1.2[round(nrow(theta.district.cov1.2)/2):nrow(theta.district.cov1.2),2,],
                          theta.district.cov1.3[round(nrow(theta.district.cov1.3)/2):nrow(theta.district.cov1.3),2,])
theta.district.cov1[,3,] <- rbind(theta.district.cov1.1[round(nrow(theta.district.cov1.1)/2):nrow(theta.district.cov1.1),3,],
                          theta.district.cov1.2[round(nrow(theta.district.cov1.2)/2):nrow(theta.district.cov1.2),3,],
                          theta.district.cov1.3[round(nrow(theta.district.cov1.3)/2):nrow(theta.district.cov1.3),3,])

sd.x.append <- c(sd.x.1[round(nrow(sd.x.1)/2):nrow(sd.x.1)],
                 sd.x.2[round(nrow(sd.x.2)/2):nrow(sd.x.2)],
                 sd.x.3[round(nrow(sd.x.3)/2):nrow(sd.x.3)])

theta.province.2.append <- rbind(theta.province.2.1[round(nrow(theta.province.2.1)/2):
                                                    nrow(theta.province.2.1), ],
                                 theta.province.2.2[round(nrow(theta.province.2.2)/2):
                                                    nrow(theta.province.2.2), ],
                                 theta.province.2.3[round(nrow(theta.province.2.3)/2):
                                                    nrow(theta.province.2.3), ])

theta.province.3.append <- rbind(theta.province.3.1[round(nrow(theta.province.3.1)/2):
                                                    nrow(theta.province.3.1), ],
                                 theta.province.3.2[round(nrow(theta.province.3.2)/2):
                                                    nrow(theta.province.3.2), ],
                                 theta.province.3.3[round(nrow(theta.province.3.3)/2):
                                                    nrow(theta.province.3.3), ])

theta.district.2.append <- rbind(theta.district.2.1[round(nrow(theta.district.2.1)/2):
                                                    nrow(theta.district.2.1), ],
                                 theta.district.2.2[round(nrow(theta.district.2.2)/2):
                                                    nrow(theta.district.2.2), ],
                                 theta.district.2.3[round(nrow(theta.district.2.3)/2):
                                                    nrow(theta.district.2.3), ])

theta.district.3.append <- rbind(theta.district.3.1[round(nrow(theta.district.3.1)/2):
                                                    nrow(theta.district.3.1), ],
                                 theta.district.3.2[round(nrow(theta.district.3.2)/2):
                                                    nrow(theta.district.3.2), ],
                                 theta.district.3.3[round(nrow(theta.district.3.3)/2):
                                                    nrow(theta.district.3.3), ])

lambda.village.2.append <- rbind(lambda.village.2.1[round(nrow(lambda.village.2.1)/2):
                                                    nrow(lambda.village.2.1), ],
                                 lambda.village.2.2[round(nrow(lambda.village.2.2)/2):
                                                    nrow(lambda.village.2.2), ],
                                 lambda.village.2.3[round(nrow(lambda.village.2.3)/2):
                                                    nrow(lambda.village.2.3), ])

lambda.village.3.append <- rbind(lambda.village.3.1[round(nrow(lambda.village.3.1)/2):
                                                    nrow(lambda.village.3.1), ],
                                 lambda.village.3.2[round(nrow(lambda.village.3.2)/2):
                                                    nrow(lambda.village.3.2), ],
                                 lambda.village.3.3[round(nrow(lambda.village.3.3)/2):
                                                    nrow(lambda.village.3.3), ])

theta.province.treat1.mean.cov <- apply(theta.province.2.append, 2, mean)
theta.province.treat1.var.cov <- apply(theta.province.2.append, 2, var)
theta.province.treat1.ci.cov <- t(apply(theta.province.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.province.treat1.ci.cov) <- c("theta.province.treat1.ci.lower.cov", "theta.province.treat1.ci.upper.cov")

theta.province.treat2.mean.cov <- apply(theta.province.3.append, 2, mean)
theta.province.treat2.var.cov <- apply(theta.province.3.append, 2, var)
theta.province.treat2.ci.cov <- t(apply(theta.province.3.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.province.treat2.ci.cov) <- c("theta.province.treat2.ci.lower.cov", "theta.province.treat2.ci.upper.cov")

theta.province.diff.mean.cov <- apply(theta.province.3.append - theta.province.2.append, 2, mean)
theta.province.diff.var.cov <- apply(theta.province.3.append - theta.province.2.append, 2, var)
theta.province.diff.ci.cov <- t(apply(theta.province.3.append - theta.province.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.province.diff.ci.cov) <- c("theta.province.diff.ci.lower.cov", "theta.province.diff.ci.upper.cov")

theta.district.treat1.mean.cov <- apply(theta.district.2.append, 2, mean)
theta.district.treat1.var.cov <- apply(theta.district.2.append, 2, var)
theta.district.treat1.ci.cov <- t(apply(theta.district.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.district.treat1.ci.cov) <- c("theta.district.treat1.ci.lower.cov", "theta.district.treat1.ci.upper.cov")

theta.district.treat2.mean.cov <- apply(theta.district.3.append, 2, mean)
theta.district.treat2.var.cov <- apply(theta.district.3.append, 2, var)
theta.district.treat2.ci.cov <- t(apply(theta.district.3.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.district.treat2.ci.cov) <- c("theta.district.treat2.ci.lower.cov", "theta.district.treat2.ci.upper.cov")

theta.district.diff.mean.cov <- apply(theta.district.3.append - theta.district.2.append, 2, mean)
theta.district.diff.var.cov <- apply(theta.district.3.append - theta.district.2.append, 2, var)
theta.district.diff.ci.cov <- t(apply(theta.district.3.append - theta.district.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(theta.district.diff.ci.cov) <- c("theta.district.diff.ci.lower.cov", "theta.district.diff.ci.upper.cov")

lambda.village.treat1.mean.cov <- apply(lambda.village.2.append, 2, mean)
lambda.village.treat1.var.cov <- apply(lambda.village.2.append, 2, var)
lambda.village.treat1.ci.cov <- t(apply(lambda.village.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(lambda.village.treat1.ci.cov) <- c("lambda.village.treat1.ci.lower.cov", "lambda.village.treat1.ci.upper.cov")

lambda.village.treat2.mean.cov <- apply(lambda.village.3.append, 2, mean)
lambda.village.treat2.var.cov <- apply(lambda.village.3.append, 2, var)
lambda.village.treat2.ci.cov <- t(apply(lambda.village.3.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(lambda.village.treat2.ci.cov) <- c("lambda.village.treat2.ci.lower.cov", "lambda.village.treat2.ci.upper.cov")

lambda.village.diff.mean.cov <- apply(lambda.village.3.append - lambda.village.2.append, 2, mean)
lambda.village.diff.var.cov <- apply(lambda.village.3.append - lambda.village.2.append, 2, var)
lambda.village.diff.ci.cov <- t(apply(lambda.village.3.append - lambda.village.2.append, 2, function(x) c(quantile(x, .025), quantile(x, .975))))
colnames(lambda.village.diff.ci.cov) <- c("lambda.village.diff.ci.lower.cov", "lambda.village.diff.ci.upper.cov")

district.lambda.mean.taliban.cov <- district.lambda.ci.lower.taliban.cov <- district.lambda.ci.upper.taliban.cov <- district.lambda.mean.isaf.cov <- district.lambda.ci.lower.isaf.cov <- district.lambda.ci.upper.isaf.cov <- rep(NA, 21)
for (dist in 1:21) {
  village.cond <- unique(data[,c("village.id", "district.id", "province.id", "M5", "Q7.1", "Q7.2", "PPLID", "DISID")])[order(unique(data[,c("village.id", "district.id", "province.id")])[,1]),]$district.id == dist
  district.lambda.mean.taliban.cov[dist] <- mean(apply(lambda.village.2.append[,village.cond], 1, mean))
  district.lambda.ci.lower.taliban.cov[dist] <- quantile(apply(lambda.village.2.append[,village.cond], 1, mean), .025)
  district.lambda.ci.upper.taliban.cov[dist] <- quantile(apply(lambda.village.2.append[,village.cond], 1, mean), .975)
  district.lambda.mean.isaf.cov[dist] <- mean(apply(lambda.village.3.append[,village.cond], 1, mean))
  district.lambda.ci.lower.isaf.cov[dist] <- quantile(apply(lambda.village.3.append[,village.cond], 1, mean), .025)
  district.lambda.ci.upper.isaf.cov[dist] <- quantile(apply(lambda.village.3.append[,village.cond], 1, mean), .975)
}


##
province.theta.mean.taliban.cov <- province.theta.ci.lower.taliban.cov <- province.theta.ci.upper.taliban.cov <- province.theta.mean.isaf.cov <- province.theta.ci.lower.isaf.cov <- province.theta.ci.upper.isaf.cov <- rep(NA, 5)
for (prov in 1:5) {
  district.cond <- unique(data[,c("district.id", "province.id", "M5", "Q7.1", "DISID")])[order(unique(data[,c("district.id", "province.id")])[,1]),]$province.id == prov
  province.theta.mean.taliban.cov[prov] <- mean(apply(theta.district.2.append[,district.cond], 1, mean))
  province.theta.ci.lower.taliban.cov[prov] <- quantile(apply(theta.district.2.append[,district.cond], 1, mean), .025)
  province.theta.ci.upper.taliban.cov[prov] <- quantile(apply(theta.district.2.append[,district.cond], 1, mean), .975)
  province.theta.mean.isaf.cov[prov] <- mean(apply(theta.district.3.append[,district.cond], 1, mean))
  province.theta.ci.lower.isaf.cov[prov] <- quantile(apply(theta.district.3.append[,district.cond], 1, mean), .025)
  province.theta.ci.upper.isaf.cov[prov] <- quantile(apply(theta.district.3.append[,district.cond], 1, mean), .975)
}

## Calculate simple difference in means estimates

taliban.district.support <- isaf.district.support <- rep(NA, length(unique(data$district.id)))
for (dist in unique(data$district.id)) {
  taliban.district.support[dist] <- mean(apply(Y[endorser[,1] == 3 & data$district.id == dist, ], 2, mean) -
                                         apply(Y[endorser[,1] == 1 & data$district.id == dist, ], 2, mean))
  ##taliban.district.support.se[dist] <- var(apply(Y[endorser[,j] == 3 & data$district.id == dist, ], 2, var) +
  ##                                       apply(Y[endorser[,j] == 1 & data$district.id == dist, ], 2, var))
  isaf.district.support[dist] <- mean(apply(Y[endorser[,1] == 2 & data$district.id == dist, ], 2, mean) -
                                      apply(Y[endorser[,1] == 1 & data$district.id == dist, ], 2, mean))
}

taliban.village.support <- isaf.village.support <- rep(NA, length(unique(data$village.id)))
for (vill in unique(data$village.id)) {
  taliban.village.support[vill] <- mean(apply(Y[endorser[,1] == 3 & data$village.id == vill, ,drop = F], 2, mean) -
                                        apply(Y[endorser[,1] == 1 & data$village.id == vill, ,drop = F], 2, mean))
  isaf.village.support[vill] <- mean(apply(Y[endorser[,1] == 2 & data$village.id == vill, , drop = F], 2, mean) -
                                     apply(Y[endorser[,1] == 1 & data$village.id == vill, , drop = F], 2, mean))
}

## put together village and district effects matrices

village.effects <- unique(data[,c("village.id", "district.id", "province.id", "M5", "Q7.1", "Q7.2", "PPLID", "DISID")])[order(unique(data[,c("village.id", "district.id", "province.id")])[,1]),]
village.effects <- cbind(village.effects, lambda.village.treat1.mean.cov, lambda.village.treat1.ci.cov, lambda.village.treat2.mean.cov, lambda.village.treat2.ci.cov, lambda.village.diff.mean.cov, lambda.village.diff.ci.cov)
names(village.effects)[4:6] <- c("province.name","district.name","village.name")
village.effects$n.village <- tapply(data$village.id, data$village.id, length)
village.effects <- village.effects[order(paste(village.effects$district.name, village.effects$village.name)),]

district.effects <- unique(data[,c("district.id", "province.id", "M5", "Q7.1", "DISID", "control")])[order(unique(data[,c("district.id", "province.id")])[,1]),]
district.effects <- cbind(district.effects, theta.district.treat1.mean.cov, theta.district.treat1.ci.cov, theta.district.treat2.mean.cov, theta.district.treat2.ci.cov, theta.district.diff.mean.cov, theta.district.diff.ci.cov)
##district.effects <- cbind(district.effects, district.lambda.mean.taliban, district.lambda.ci.lower.taliban, district.lambda.ci.upper.taliban, district.lambda.mean.isaf, district.lambda.ci.lower.isaf, district.lambda.ci.upper.isaf, district.lambda.mean.taliban.cov, district.lambda.ci.lower.taliban.cov, district.lambda.ci.upper.taliban.cov, district.lambda.mean.isaf.cov, district.lambda.ci.lower.isaf.cov, district.lambda.ci.upper.isaf.cov)
district.effects$n.district <- tapply(data$province.id, data$district.id, length)
names(district.effects)[3:4] <- c("province.name","district.name")
district.effects <- district.effects[order(paste(district.effects$province.name, district.effects$district.name)),]

province.effects <- unique(data[,c("province.id", "M5")])[order(unique(data[,c("province.id")])),]
province.effects <- cbind(province.effects, theta.province.treat1.mean.cov, theta.province.treat1.ci.cov, theta.province.treat2.mean.cov, theta.province.treat2.ci.cov, theta.province.diff.mean.cov, theta.province.diff.ci.cov)
province.effects <- cbind(province.effects, province.theta.mean.taliban.cov, province.theta.ci.lower.taliban.cov, province.theta.ci.upper.taliban.cov, province.theta.mean.isaf.cov, province.theta.ci.lower.isaf.cov, province.theta.ci.upper.isaf.cov)
province.effects$n.province <- tapply(data$province.id, data$province.id, length)
names(province.effects)[2] <- c("province.name")
province.effects <- province.effects[order(paste(province.effects$province.name, province.effects$province.name)),]


## By district plot

district.counts <- tapply(district.effects$district.id, district.effects$province.id, length)



### Now with covariate estimates

## By district plot

district.counts <- tapply(district.effects$district.id, district.effects$province.id, length)

## with diff in means est

if(!exists("no.covariate.model")) {
  add.text <- "."
} else {
  add.text <- ""
}

if (!exists("no.covariate.model")) {
  if(!exists("basic.model")) {
    tb.indiv.treat1 <- round(cbind(apply(lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,2,], 2, sd)),3)
    tb.village.treat1 <- round(cbind(apply(theta.village.cov1[,2,], 2, mean), apply(theta.village.cov1[,2,], 2, sd)),3)
    tb.district.treat1 <- round(cbind(apply(theta.district.cov1[,2,], 2, mean), apply(theta.district.cov1[,2,], 2, sd)),3)
    tb.indiv.treat2 <- round(cbind(apply(lambda.cov1[,3,], 2, mean), apply(lambda.cov1[,3,], 2, sd)),3)
    tb.village.treat2 <- round(cbind(apply(theta.village.cov1[,3,], 2, mean), apply(theta.village.cov1[,3,], 2, sd)),3)
    tb.district.treat2 <- round(cbind(apply(theta.district.cov1[,3,], 2, mean), apply(theta.district.cov1[,3,], 2, sd)),3)
    
    rownames(tb.indiv.treat1) <- paste("indv.isaf.",colnames(indivcov),sep="")
    rownames(tb.indiv.treat2) <- paste("indv.talb.",colnames(indivcov),sep="")
    rownames(tb.village.treat1) <-paste("vill.isaf.", colnames(villcov), sep ="")
    rownames(tb.village.treat2) <- paste("vill.talb.", colnames(villcov), sep ="")
    rownames(tb.district.treat1) <-paste("dist.isaf.",colnames(distcov), sep = "")
    rownames(tb.district.treat2) <- paste("dist.talb.",colnames(distcov), sep = "")

  
    tb <- rbind(tb.indiv.treat1, tb.indiv.treat2, tb.village.treat1, tb.village.treat2, tb.district.treat1, tb.district.treat2)
    colnames(tb)<-c("est","se")
    write.csv(tb, file = paste("results/estimates.",file.name, ".csv",sep=""))
    
  } else {
    tb.indiv.treat1 <- round(cbind(apply(lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,2,], 2, sd)),3)
    tb.indiv.treat2 <- round(cbind(apply(lambda.cov1[,3,], 2, mean), apply(lambda.cov1[,3,], 2, sd)),3)
    
    rownames(tb.indiv.treat1) <- paste("indv.isaf.",colnames(indivcov),sep="")
    rownames(tb.indiv.treat2) <- paste("indv.talb.",colnames(indivcov),sep="")
  
    tb <- rbind(tb.indiv.treat1, tb.indiv.treat2)
    colnames(tb)<-c("est","se")
    write.csv(tb, file = paste("results/estimates.",file.name, ".csv",sep=""))
  }
}

## now write them with 95% CI's and not rounded
if (!exists("no.covariate.model")) {
  if(!exists("basic.model")) {
    tb.indiv.treat1 <- cbind(apply(lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,2,], 2, sd),
                             apply(lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.village.treat1 <- cbind(apply(theta.village.cov1[,2,], 2, mean), apply(theta.village.cov1[,2,], 2, sd),
                               apply(theta.village.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.district.treat1 <- cbind(apply(theta.district.cov1[,2,], 2, mean), apply(theta.district.cov1[,2,], 2, sd),
                                apply(theta.district.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.indiv.treat2 <- cbind(apply(lambda.cov1[,3,], 2, mean), apply(lambda.cov1[,3,], 2, sd),
                             apply(lambda.cov1[,3,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,], 2, function(x) quantile(x, .975)))
    tb.village.treat2 <- cbind(apply(theta.village.cov1[,3,], 2, mean), apply(theta.village.cov1[,3,], 2, sd),
                               apply(theta.village.cov1[,3,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,3,], 2, function(x) quantile(x, .975)))
    tb.district.treat2 <- cbind(apply(theta.district.cov1[,3,], 2, mean), apply(theta.district.cov1[,3,], 2, sd),
                                apply(theta.district.cov1[,3,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,3,], 2, function(x) quantile(x, .975)))

    tb.indiv.diff <- cbind(apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, sd),
                             apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.village.diff <- cbind(apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, mean), apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, sd),
                               apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.district.diff <- cbind(apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, mean), apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, sd),
                                apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, function(x) quantile(x, .975)))
    
    rownames(tb.indiv.treat1) <- paste("indv.isaf.",colnames(indivcov),sep="")
    rownames(tb.indiv.treat2) <- paste("indv.talb.",colnames(indivcov),sep="")
    rownames(tb.indiv.diff) <- paste("indv.diff.",colnames(indivcov),sep="")
    rownames(tb.village.treat1) <-paste("vill.isaf.", colnames(villcov), sep ="")
    rownames(tb.village.treat2) <- paste("vill.talb.", colnames(villcov), sep ="")
    rownames(tb.village.diff) <- paste("vill.diff.", colnames(villcov), sep ="")
    rownames(tb.district.treat1) <-paste("dist.isaf.",colnames(distcov), sep = "")
    rownames(tb.district.treat2) <- paste("dist.talb.",colnames(distcov), sep = "")
    rownames(tb.district.diff) <- paste("dist.diff.",colnames(distcov), sep = "")
  
    tb <- rbind(tb.indiv.treat1, tb.indiv.treat2, tb.indiv.diff, tb.village.treat1, tb.village.treat2, tb.village.diff, tb.district.treat1, tb.district.treat2, tb.district.diff)
    colnames(tb)<-c("est","se","lower.ci","upper.ci")
    write.csv(tb, file = paste("results/estimates.ci.",file.name, ".csv",sep=""))

    if(!exists("harm.only")) {
      tb.harmtaliban.treat1 <- cbind(mean(lambda.cov1[,2,1] + lambda.cov1[,2,5]),
                                     sd(lambda.cov1[,2,1] + lambda.cov1[,2,5]),
                                     quantile(lambda.cov1[,2,1] + lambda.cov1[,2,5], .025),
                                     quantile(lambda.cov1[,2,1] + lambda.cov1[,2,5], .975))
      tb.harmISAF.treat1 <- cbind(mean(lambda.cov1[,2,3] + lambda.cov1[,2,7]),
                                  sd(lambda.cov1[,2,3] + lambda.cov1[,2,7]),
                                  quantile(lambda.cov1[,2,3] + lambda.cov1[,2,7], .025),
                                  quantile(lambda.cov1[,2,3] + lambda.cov1[,2,7], .975))
      
      tb.harmtaliban.treat2 <- cbind(mean(lambda.cov1[,3,1] + lambda.cov1[,3,5]),
                                     sd(lambda.cov1[,3,1] + lambda.cov1[,3,5]),
                                     quantile(lambda.cov1[,3,1] + lambda.cov1[,3,5], .025),
                                     quantile(lambda.cov1[,3,1] + lambda.cov1[,3,5], .975))
      tb.harmISAF.treat2 <- cbind(mean(lambda.cov1[,3,3] + lambda.cov1[,3,7]),
                                  sd(lambda.cov1[,3,3] + lambda.cov1[,3,7]),
                                  quantile(lambda.cov1[,3,3] + lambda.cov1[,3,7], .025),
                                  quantile(lambda.cov1[,3,3] + lambda.cov1[,3,7], .975))
      
      tb.harmtaliban.diff <- cbind(mean(lambda.cov1[,3,1] + lambda.cov1[,3,5] - (lambda.cov1[,2,1] + lambda.cov1[,2,5])),
                                   sd(lambda.cov1[,3,1] + lambda.cov1[,3,5] - (lambda.cov1[,2,1] + lambda.cov1[,2,5])),
                                   quantile(lambda.cov1[,3,1] + lambda.cov1[,3,5] - (lambda.cov1[,2,1] + lambda.cov1[,2,5]), .025),
                                   quantile(lambda.cov1[,3,1] + lambda.cov1[,3,5] - (lambda.cov1[,2,1] + lambda.cov1[,2,5]), .975))
      tb.harmISAF.diff <- cbind(mean(lambda.cov1[,3,3] + lambda.cov1[,3,7] - (lambda.cov1[,2,3] + lambda.cov1[,2,7])),
                                sd(lambda.cov1[,3,3] + lambda.cov1[,3,7] - (lambda.cov1[,2,3] + lambda.cov1[,2,7])),
                                quantile(lambda.cov1[,3,3] + lambda.cov1[,3,7] - (lambda.cov1[,2,3] + lambda.cov1[,2,7]), .025),
                                quantile(lambda.cov1[,3,3] + lambda.cov1[,3,7] - (lambda.cov1[,2,3] + lambda.cov1[,2,7]), .975))
      
      tb <- rbind(tb.harmtaliban.treat1, tb.harmtaliban.treat2, tb.harmtaliban.diff, tb.harmISAF.treat1, tb.harmISAF.treat2, tb.harmISAF.diff)
      rownames(tb) <- c("indv.harmtalib.isaf", "indv.harmtalib.taliban", "indv.harmtalib.diff","indv.harmISAF.isaf", "indv.harmISAF.taliban","indv.harmISAF.diff")
      colnames(tb)<-c("est","se","lower.ci","upper.ci")
      write.csv(tb, file = paste("results/estimatesHarmApproach.ci.",file.name, ".csv",sep=""))
    }
    
  } else {
    tb.indiv.treat1 <- cbind(apply(lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,2,], 2, sd),
                             apply(lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
    tb.indiv.treat2 <- cbind(apply(lambda.cov1[,3,], 2, mean), apply(lambda.cov1[,3,], 2, sd),
                             apply(lambda.cov1[,3,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,], 2, function(x) quantile(x, .975)))
    tb.indiv.diff <- cbind(apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, sd),
                             apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
  
    rownames(tb.indiv.treat1) <- paste("indv.isaf.",colnames(indivcov),sep="")
    rownames(tb.indiv.treat2) <- paste("indv.talb.",colnames(indivcov),sep="")
    rownames(tb.indiv.diff) <- paste("indv.diff.",colnames(indivcov),sep="")
    
    tb <- rbind(tb.indiv.treat1, tb.indiv.treat2)
    colnames(tb)<-c("est","se","lower.ci","upper.ci")
    write.csv(tb, file = paste("results/estimates.ci.",file.name, ".csv",sep=""))
  }
}


if(exists("pro.taliban.model") == TRUE | exists("pro.taliban.approach.model") | exists("prov.id") == TRUE) {  
  
  tb.indiv.treat1 <- cbind(apply(lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,2,], 2, sd),
                           apply(lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
  tb.village.treat1 <- cbind(apply(theta.village.cov1[,2,], 2, mean), apply(theta.village.cov1[,2,], 2, sd),
                             apply(theta.village.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,2,], 2, function(x) quantile(x, .975)))
  tb.district.treat1 <- cbind(apply(theta.district.cov1[,2,], 2, mean), apply(theta.district.cov1[,2,], 2, sd),
                              apply(theta.district.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,2,], 2, function(x) quantile(x, .975)))
  tb.indiv.treat2 <- cbind(apply(lambda.cov1[,3,], 2, mean), apply(lambda.cov1[,3,], 2, sd),
                           apply(lambda.cov1[,3,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,], 2, function(x) quantile(x, .975)))
  tb.village.treat2 <- cbind(apply(theta.village.cov1[,3,], 2, mean), apply(theta.village.cov1[,3,], 2, sd),
                             apply(theta.village.cov1[,3,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,3,], 2, function(x) quantile(x, .975)))
  tb.district.treat2 <- cbind(apply(theta.district.cov1[,3,], 2, mean), apply(theta.district.cov1[,3,], 2, sd),
                              apply(theta.district.cov1[,3,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,3,], 2, function(x) quantile(x, .975)))
  
  tb.indiv.diff <- cbind(apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, mean), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, sd),
                         apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .025)), apply(lambda.cov1[,3,] - lambda.cov1[,2,], 2, function(x) quantile(x, .975)))
  tb.village.diff <- cbind(apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, mean), apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, sd),
                           apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.village.cov1[,3,] - theta.village.cov1[,2,], 2, function(x) quantile(x, .975)))
  tb.district.diff <- cbind(apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, mean), apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, sd),
                            apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, function(x) quantile(x, .025)), apply(theta.district.cov1[,3,] - theta.district.cov1[,2,], 2, function(x) quantile(x, .975)))

  if(exists("prov.id")) {
    add.pro.taliban <- 1
  }else{
    add.pro.taliban <- 0
  }
  
  tb.violencetaliban.treat1 <- cbind(mean(lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban]),
                           sd(lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban]),
                           quantile(lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban], .025),
                           quantile(lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban], .975))
  tb.violenceISAF.treat1 <- cbind(mean(lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban]),
                           sd(lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban]),
                           quantile(lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban], .025),
                           quantile(lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban], .975))

  tb.violencetaliban.treat2 <- cbind(mean(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban]),
                           sd(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban]),
                           quantile(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban], .025),
                           quantile(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban], .975))
  tb.violenceISAF.treat2 <- cbind(mean(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban]),
                           sd(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban]),
                           quantile(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban], .025),
                           quantile(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban], .975))

  tb.violencetaliban.diff <- cbind(mean(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban] - (lambda.cov1[,2,1] + lambda.cov1[,2,15] + lambda.cov1[,2,16+add.pro.taliban])),
                           sd(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban] - (lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban])),
                           quantile(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban] - (lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban]), .025),
                           quantile(lambda.cov1[,3,1] + lambda.cov1[,3,16+add.pro.taliban] - (lambda.cov1[,2,1] + lambda.cov1[,2,16+add.pro.taliban]), .975))
  tb.violenceISAF.diff <- cbind(mean(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban] - (lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban])),
                           sd(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban] - (lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban])),
                           quantile(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban] - (lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban]), .025),
                           quantile(lambda.cov1[,3,3] + lambda.cov1[,3,17+add.pro.taliban] - (lambda.cov1[,2,3] + lambda.cov1[,2,17+add.pro.taliban]), .975))

  if(exists("pro.taliban.approach.model") | exists("prov.id")) {
    tb.approachtaliban.treat1 <- cbind(mean(lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban]),
                                       sd(lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban]),
                                       quantile(lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban], .025),
                                       quantile(lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban], .975))
    tb.approachISAF.treat1 <- cbind(mean(lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban]),
                                    sd(lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban]),
                                    quantile(lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban], .025),
                                    quantile(lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban], .975))

    tb.approachtaliban.treat2 <- cbind(mean(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban]),
                                       sd(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban]),
                                       quantile(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban], .025),
                                       quantile(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban], .975))
    tb.approachISAF.treat2 <- cbind(mean(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban]),
                                    sd(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban]),
                                    quantile(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban], .025),
                                    quantile(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban], .975))
    
    tb.approachtaliban.diff <- cbind(mean(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban] - (lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban])),
                                     sd(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban] - (lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban])),
                                     quantile(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban] - (lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban]), .025),
                                     quantile(lambda.cov1[,3,5] + lambda.cov1[,3,18+add.pro.taliban] - (lambda.cov1[,2,5] + lambda.cov1[,2,18+add.pro.taliban]), .975))
    tb.approachISAF.diff <- cbind(mean(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban] - (lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban])),
                                  sd(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban] - (lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban])),
                                  quantile(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban] - (lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban]), .025),
                                  quantile(lambda.cov1[,3,7] + lambda.cov1[,3,19+add.pro.taliban] - (lambda.cov1[,2,7] + lambda.cov1[,2,19+add.pro.taliban]), .975))
  }
  
  rownames(tb.indiv.treat1) <- paste("indv.isaf.",colnames(indivcov),sep="")
  rownames(tb.indiv.treat2) <- paste("indv.talb.",colnames(indivcov),sep="")
  rownames(tb.indiv.diff) <- paste("indv.diff.",colnames(indivcov),sep="")
  rownames(tb.village.treat1) <-paste("vill.isaf.", colnames(villcov), sep ="")
  rownames(tb.village.treat2) <- paste("vill.talb.", colnames(villcov), sep ="")
  rownames(tb.village.diff) <- paste("vill.diff.", colnames(villcov), sep ="")
  rownames(tb.district.treat1) <-paste("dist.isaf.",colnames(distcov), sep = "")
  rownames(tb.district.treat2) <- paste("dist.talb.",colnames(distcov), sep = "")
  rownames(tb.district.diff) <- paste("dist.diff.",colnames(distcov), sep = "")

  tb <- rbind(tb.indiv.treat1, tb.indiv.treat2, tb.indiv.diff, tb.village.treat1, tb.village.treat2, tb.village.diff, tb.district.treat1, tb.district.treat2, tb.district.diff,tb.violencetaliban.treat1, tb.violencetaliban.treat2, tb.violencetaliban.diff, tb.violenceISAF.treat1, tb.violenceISAF.treat2, tb.violenceISAF.diff)
  if(exists("pro.taliban.model") | exists("prov.id"))
    rownames(tb)[(length(rownames(tb))-5):length(rownames(tb))] <- c("indv.violencetalib.isaf", "indv.violencetalib.taliban", "indv.violencetalib.diff","indv.violenceISAF.isaf", "indv.violenceISAF.taliban","indv.violenceISAF.diff")
  if(exists("pro.taliban.approach.model") | exists("prov.id")) {
    tb <- rbind(tb,tb.approachtaliban.treat1, tb.approachtaliban.treat2, tb.approachtaliban.diff, tb.approachISAF.treat1, tb.approachISAF.treat2, tb.approachISAF.diff)
    rownames(tb)[(length(rownames(tb))-11):length(rownames(tb))] <-  c("indv.violencetalib.isaf", "indv.violencetalib.taliban", "indv.violencetalib.diff","indv.violenceISAF.isaf", "indv.violenceISAF.taliban","indv.violenceISAF.diff", "indv.approachtalib.isaf", "indv.approachtalib.taliban", "indv.approachtalib.diff","indv.approachISAF.isaf", "indv.approachISAF.taliban","indv.approachISAF.diff")
  }
  colnames(tb)<-c("est","se","lower.ci","upper.ci")
  write.csv(tb, file = paste("results/estimatesInteraction.ci.",file.name, ".csv",sep=""))
  
}

if(model.iter == 1) {

## make district-level difference plot

  village.effects <- unique(data[,c("village.id", "district.id", "province.id", "M5", "Q7.1", "Q7.2", "PPLID", "DISID", "violence.count.taliban.5km.std", "violence.count.ISAF.5km.std","violence.count.taliban.5km", "violence.count.ISAF.5km","altitude.std", "population.std")])[order(unique(data[,c("village.id", "district.id", "province.id")])[,1]),]
village.effects <- cbind(village.effects, lambda.village.treat1.mean.cov, lambda.village.treat1.var.cov, lambda.village.treat1.ci.cov, lambda.village.treat2.mean.cov, lambda.village.treat2.var.cov, lambda.village.treat2.ci.cov, lambda.village.diff.mean.cov, lambda.village.diff.var.cov, lambda.village.diff.ci.cov)
names(village.effects)[4:6] <- c("province.name","district.name","village.name")
village.effects$n.village <- tapply(data$village.id, data$village.id, length)
village.effects <- village.effects[order(paste(village.effects$district.name, village.effects$village.name)),]


district.effects <- unique(data[,c("district.id", "province.id", "M5", "Q7.1", "DISID", "sharia.courts", "cerp.money.std", "opium.ha.std", "cdc.projects", "roads.km.std", "pakistan", "control")])[order(unique(data[,c("district.id", "province.id")])[,1]),]
district.effects <- cbind(district.effects, theta.district.treat1.mean.cov, theta.district.treat1.var.cov, theta.district.treat1.ci.cov, theta.district.treat2.mean.cov, theta.district.treat2.var.cov, theta.district.treat2.ci.cov, theta.district.diff.mean.cov, theta.district.diff.var.cov, theta.district.diff.ci.cov)
##district.effects <- cbind(district.effects, district.lambda.mean.taliban, district.lambda.ci.lower.taliban, district.lambda.ci.upper.taliban, district.lambda.mean.isaf, district.lambda.ci.lower.isaf, district.lambda.ci.upper.isaf, district.lambda.mean.taliban.cov, district.lambda.ci.lower.taliban.cov, district.lambda.ci.upper.taliban.cov, district.lambda.mean.isaf.cov, district.lambda.ci.lower.isaf.cov, district.lambda.ci.upper.isaf.cov)
district.effects$n.district <- tapply(data$province.id, data$district.id, length)
names(district.effects)[3:4] <- c("province.name","district.name")
district.effects <- district.effects[order(paste(district.effects$province.name, district.effects$district.name)),]


dist.violence.ISAF <- aggregate(village.effects$violence.count.ISAF.5km, list(village.effects$district.id), FUN = sum)
names(dist.violence.ISAF) <- c("district.id", "violence.count.ISAF")
dist.violence.taliban <- aggregate(village.effects$violence.count.taliban.5km, list(village.effects$district.id), FUN = sum)
names(dist.violence.taliban) <- c("district.id", "violence.count.taliban")

district.effects <- merge(district.effects, dist.violence.ISAF, by = "district.id", sort = F)
district.effects <- merge(district.effects, dist.violence.taliban, by = "district.id", sort = F)

### Control

##control.taliban == 1
##control.contested == 1
##cerp,cdc,control,opium,sharia,pakistan

var.list <- list()
var.list[[3]] <- ifelse(district.effects$control == 0, 1, ifelse(district.effects$control == 1, 0, 5))

## Opium

var.list[[4]] <- district.effects$opium.ha.std >= 0

## CERP money

var.list[[1]] <- district.effects$cerp.money.std >= 0

## pakistan

##pakistan 1 / 0

var.list[[6]] <- district.effects$pakistan

## CERP money

var.list[[2]] <- district.effects$cdc.projects >= 3

## Sharia

var.list[[5]] <- district.effects$sharia.courts > 0

var.names <- c("cerp","cdc","control","opium", "sharia", "pakistan")

dist.ISAF.var.est <- dist.ISAF.var.lower.ci <- dist.ISAF.var.upper.ci<- 
  dist.ISAF.novar.est <- dist.ISAF.novar.lower.ci <- dist.ISAF.novar.upper.ci<- 
  dist.ISAF.diffvar.est <- dist.ISAF.diffvar.lower.ci <- dist.ISAF.diffvar.upper.ci<- 
  dist.taliban.var.est <- dist.taliban.var.lower.ci <- dist.taliban.var.upper.ci<- 
  dist.taliban.novar.est <- dist.taliban.novar.lower.ci <- dist.taliban.novar.upper.ci<- 
  dist.taliban.diffvar.est <- dist.taliban.diffvar.lower.ci <- dist.taliban.diffvar.upper.ci<- 
  dist.diff.var.est <- dist.diff.var.lower.ci <- dist.diff.var.upper.ci<- 
  dist.diff.novar.est <- dist.diff.novar.lower.ci <- dist.diff.novar.upper.ci<- 
  dist.diff.diffvar.est <- dist.diff.diffvar.lower.ci <- dist.diff.diffvar.upper.ci <- list()
      
for (i in 1:length(var.list)) {
  
  dist.ISAF.var.est[[var.names[i]]] <- mean(apply(theta.district.2.append[,var.list[[i]]==1], 1, mean))
  dist.ISAF.var.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==1], 1, mean), .025)
  dist.ISAF.var.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==1], 1, mean), .975)
  
  dist.ISAF.novar.est[[var.names[i]]] <- mean(apply(theta.district.2.append[,var.list[[i]]==0], 1, mean))
  dist.ISAF.novar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==0], 1, mean), .025)
  dist.ISAF.novar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==0], 1, mean), .975)
  
  dist.ISAF.diffvar.est[[var.names[i]]] <- mean(apply(theta.district.2.append[,var.list[[i]]==1], 1, mean) - apply(theta.district.2.append[,var.list[[i]]==0], 1, mean))
  dist.ISAF.diffvar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==1],1,mean) - apply(theta.district.2.append[,var.list[[i]]==0], 1, mean), .025)
  dist.ISAF.diffvar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.2.append[,var.list[[i]]==1],1,mean) - apply(theta.district.2.append[,var.list[[i]]==0], 1, mean), .975)
  
  dist.taliban.var.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean))
  dist.taliban.var.lower.ci[[var.names[i]]] <-  quantile(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean), .025)
  dist.taliban.var.upper.ci[[var.names[i]]] <-  quantile(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean), .975)
  
  dist.taliban.novar.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==0], 1, mean))
  dist.taliban.novar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==0], 1, mean), .025)
  dist.taliban.novar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==0], 1, mean), .975)
  
  dist.taliban.diffvar.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean) - apply(theta.district.3.append[,var.list[[i]]==0], 1, mean))
  dist.taliban.diffvar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean) - apply(theta.district.3.append[,var.list[[i]]==0], 1, mean), .025)
  dist.taliban.diffvar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==1], 1, mean) - apply(theta.district.3.append[,var.list[[i]]==0], 1, mean), .975)
  
  dist.diff.var.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean))
  dist.diff.var.lower.ci[[var.names[i]]] <-  quantile(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean), .025)
  dist.diff.var.upper.ci[[var.names[i]]] <-  quantile(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean), .975)
  
  dist.diff.novar.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean))
  dist.diff.novar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean), .025)
  dist.diff.novar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean), .975)
  
  dist.diff.diffvar.est[[var.names[i]]] <- mean(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean) -
                                  apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean))
  dist.diff.diffvar.lower.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean) -
                                           apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean), .025)
  dist.diff.diffvar.upper.ci[[var.names[i]]] <- quantile(apply(theta.district.3.append[,var.list[[i]]==1] - theta.district.2.append[,var.list[[i]]==1], 1, mean) -
                                           apply(theta.district.3.append[,var.list[[i]]==0] - theta.district.2.append[,var.list[[i]]==0], 1, mean), .975)
}

pdf("figs/Fig12-DistrictDiffMeans.pdf", height = 5, width = 6.5)
par(mfcol = c(3,1), las = 1, mar=c(0,1.5,0,1), oma=c(0,4,0,0), cex = 0.6)

x.coords <- c(1.33, 1.66, 1.99,   2.66, 2.99, 3.33,   3.99, 4.33, 4.66,   5.33, 5.66, 5.99,  6.66, 6.99, 7.33,   7.99, 8.33, 8.66)

### ISAF Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1.23, 8.76), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Support for Taliban", outer = TRUE, side = 2, cex = 0.65, at = 0.85, line = 2.5)

points(x.coords, as.vector(rbind(do.call(c, dist.taliban.var.est), do.call(c, dist.taliban.novar.est),
                                 do.call(c, dist.taliban.diffvar.est))),
       pch = rep(c(19,15,2), 2), cex = 1.5)
segments(x.coords, as.vector(rbind(do.call(c, dist.taliban.var.lower.ci), do.call(c, dist.taliban.novar.lower.ci),
                                   do.call(c, dist.taliban.diffvar.lower.ci))),
         x.coords, as.vector(rbind(do.call(c, dist.taliban.var.upper.ci), do.call(c, dist.taliban.novar.upper.ci),
                                   do.call(c, dist.taliban.diffvar.upper.ci))))


text(x.coords[1], -1.2, "High", cex = .75)
text(x.coords[2], -1.2, "Low", cex = .75)
text(x.coords[3], -.6, "Difference\nHigh - Low", cex = .75)

text(x.coords[4], -1.2, "0 - 2", cex = .75)
text(x.coords[5], -1.2, "3 - 6", cex = .75)
text(x.coords[6], -0.6, "Difference\nMany - Few", cex = .75)

text(x.coords[7], -0.9, "Taliban", cex = .75)
text(x.coords[8], -1.2, "Contested", cex = .75)
text(x.coords[9], -.6, "Difference\nTaliban-\nContested", cex = .75)

text(x.coords[10], -1.2, "High", cex = .75)
text(x.coords[11], -1.2, "Low", cex = .75)
text(x.coords[12], -0.8, "Difference\nHigh - Low", cex = .75)
text(x.coords[13], -1.2, "Yes", cex = .75)
text(x.coords[14], -1.2, "No", cex = .75)
text(x.coords[15], -.8, "Difference\nYes - No", cex = .75)

text(x.coords[16], -1.2, "Yes", cex = .75)
text(x.coords[17], -1.2, "No", cex = .75)
text(x.coords[18], -.8, "Difference\nYes - No", cex = .75)

mtext("CERP Spending", outer = TRUE, cex = 0.6, at = 0.115 + .155*0, line = -1.5)
mtext("CDC Projects", outer = TRUE, cex = 0.6, at = 0.115 + .155*1, line = -1.5)
mtext("Territorial Control", outer = TRUE, cex = 0.6, at = 0.115 + .155*2, line = -1.5)
mtext("Opium Production", outer = TRUE, cex = 0.6, at = 0.115 + .155*3, line = -1.5)
mtext("Sharia Courts", outer = TRUE, cex = 0.6, at = 0.115 + .155*4, line = -1.5)
mtext("Borders Pakistan", outer = TRUE, cex = 0.6, at = 0.115 + .155*5, line = -1.5)

### Taliban Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1.23, 8.76), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Support for ISAF", outer = TRUE, side = 2, cex = 0.65, at = 0.51, line = 2.5)


points(x.coords, as.vector(rbind(do.call(c, dist.ISAF.var.est), do.call(c, dist.ISAF.novar.est),
                                 do.call(c, dist.ISAF.diffvar.est))),
       pch = rep(c(19,15,2), 2), cex = 1.5)
segments(x.coords, as.vector(rbind(do.call(c, dist.ISAF.var.lower.ci), do.call(c, dist.ISAF.novar.lower.ci),
                                   do.call(c, dist.ISAF.diffvar.lower.ci))),
         x.coords, as.vector(rbind(do.call(c, dist.ISAF.var.upper.ci), do.call(c, dist.ISAF.novar.upper.ci),
                                   do.call(c, dist.ISAF.diffvar.upper.ci))))

### Difference Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1.23, 8.76), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Difference in Support\n Taliban - ISAF", outer = TRUE, side = 2, padj = .6, cex = 0.65, at = 0.16, line = 2.5)

points(x.coords, as.vector(rbind(do.call(c, dist.diff.var.est), do.call(c, dist.diff.novar.est),
                                 do.call(c, dist.diff.diffvar.est))),
       pch = rep(c(19,15,2), 2), cex = 1.5)
segments(x.coords, as.vector(rbind(do.call(c, dist.diff.var.lower.ci), do.call(c, dist.diff.novar.lower.ci),
                                   do.call(c, dist.diff.diffvar.lower.ci))),
         x.coords, as.vector(rbind(do.call(c, dist.diff.var.upper.ci), do.call(c, dist.diff.novar.upper.ci),
                                   do.call(c, dist.diff.diffvar.upper.ci))))

dev.off()
  
## make harm descriptive plot
  
cov <- as.data.frame(indivcov)

cov$violent.exp.taliban <- 1
cov$violent.exp.ISAF <- 0
talibanharm.ISAF <- apply(lambda.cov1[ , 2, ] %*% t(cov) + lambda.village.2.append[, village], 1, mean)

cov$violent.exp.taliban <- 0
cov$violent.exp.ISAF <- 1
ISAFharm.ISAF <- apply(lambda.cov1[ , 2, ] %*% t(cov) + lambda.village.2.append[, village], 1, mean)

cov$violent.exp.taliban <- 0
cov$violent.exp.ISAF <- 0
noharm.ISAF <- apply(lambda.cov1[ , 2, ] %*% t(cov) + lambda.village.2.append[, village], 1, mean)

cov$violent.exp.taliban <- 1
cov$violent.exp.ISAF <- 1
bothharm.ISAF <- apply(lambda.cov1[ , 2, ] %*% t(cov) + lambda.village.2.append[, village], 1, mean)


cov$violent.exp.taliban <- 1
cov$violent.exp.ISAF <- 0
talibanharm.taliban <- apply(lambda.cov1[ , 3, ] %*% t(cov) + lambda.village.3.append[, village], 1, mean)

cov$violent.exp.taliban <- 0
cov$violent.exp.ISAF <- 1
ISAFharm.taliban <- apply(lambda.cov1[ , 3, ] %*% t(cov) + lambda.village.3.append[, village], 1, mean)

cov$violent.exp.taliban <- 0
cov$violent.exp.ISAF <- 0
noharm.taliban <- apply(lambda.cov1[ , 3, ] %*% t(cov) + lambda.village.3.append[, village], 1, mean)

cov$violent.exp.taliban <- 1
cov$violent.exp.ISAF <- 1
bothharm.taliban <- apply(lambda.cov1[ , 3, ] %*% t(cov) + lambda.village.3.append[, village], 1, mean)

summ <- function(x) c(mean(x), quantile(x, .025), quantile(x, .975))

pdf("figs/Fig5-CoefsHarmDescriptive.pdf", height = 5/3, width = 5)

par(mfcol = c(1,2), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1,2,3,4)

### ISAF Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.1, 1), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Support for ISAF", outer = FALSE, side = 2, cex = 0.65, at = 0, line = 3, padj = 0)

coefsPlot <- rbind(summ(noharm.ISAF), summ(talibanharm.ISAF), summ(ISAFharm.ISAF), summ(bothharm.ISAF) )

points(x.coords, coefsPlot[,1], pch = rep(c(19), 4), cex = 1.5)
segments(x.coords, coefsPlot[,2], x.coords, coefsPlot[,3])

text(x.coords[1] , 1.1, c("None"), pos = 1)
text(x.coords[2] , 1.1, c("Taliban\nOnly"), pos = 1)
text(x.coords[3] , 1.1, c("ISAF\nOnly"), pos = 1)
text(x.coords[4] , 1.1, c("Both"), pos = 1)

mtext("Victimization", outer = TRUE, cex = 0.65, at = .225, line = -0.5, padj = -1.4)
##mtext("by ISAF", outer = TRUE, cex = 0.65, at = .43, line = -0.5)
##mtext("by Taliban", outer = TRUE, cex = 0.65, at = .54, line = -0.5)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.1, 1), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Support for the Taliban", outer = FALSE, side = 2, cex = 0.65, at = 0, line = 3, padj = 0)

coefsPlot <- rbind(summ(noharm.taliban), summ(talibanharm.taliban), summ(ISAFharm.taliban), summ(bothharm.taliban) )

points(x.coords, coefsPlot[,1], pch = rep(c(19), 4), cex = 1.5)
segments(x.coords, coefsPlot[,2], x.coords, coefsPlot[,3])

text(x.coords[1] , 1.1, c("None"), pos = 1)
text(x.coords[2] , 1.1, c("Taliban\nOnly"), pos = 1)
text(x.coords[3] , 1.1, c("ISAF\nOnly"), pos = 1)
text(x.coords[4] , 1.1, c("Both"), pos = 1)

mtext("Victimization", outer = TRUE, cex = 0.65, at = .73, line = -0.5, padj = -1.4)

##text(x.coords[4] , -.5, "Manteqa\n(area)", cex = 1)

##mtext("Victimization", outer = TRUE, cex = 0.65, at = 0.16, line = -0.5, padj = -2)
##mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.1, line = -0.5)
##mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.21, line = -0.5)

dev.off()
}

}




## first plot the key violence vars for the full model

library(xtable)

####
### new version of harm plot

indivCoefs <- read.csv("results/estimates.ci..individual.fit.csv")
areaCoefs <- read.csv("results/estimates.ci..area.fit.csv")

tabISAF <- cbind(indivCoefs[c(1:16, 49:52, 61:68),2:3], areaCoefs[c(1:16, 49:52, 61:68),2:3])
tabTaliban <- cbind(indivCoefs[c(17:32, 53:56, 69:76),2:3], areaCoefs[c(17:32, 53:56, 69:76),2:3])

rownames(tabISAF) <- rownames(tabTaliban) <- c("Harm from Taliban violence", "Harm from Taliban violence is NA",
                   "Harm from ISAF violence", "Harm from ISAF violence is NA",
                   "Approach by Taliban after Harm", "Approach by Taliban after Harm is NA",
                   "Approach by ISAF after Harm", "Approach by ISAF after Harm is NA",
                   "ISAF encounter frequency", "Years of education",
                   "Age (tens)", "Income (Afghanis)", "Income is NA",
                   "Schooled in madrassa", "Pro-Taliban tribe",
                   "Pro-Taliban tribe is NA", "Altitude (km)",
                   "Population", "ISAF-initiated violent events (within 5km)",
                   "Taliban-initiated violent events (within 5km)",
                   "Sha'ria courts", "CERP project spending",
                   "Opium cultivation (ha.)", "CDC project count",
                   "Road length (km)", "Pakistan border",
                   "Government territorial control","Contested territorial control")
                   
sink("figs/afghanRegResults.txt")
print("ISAF")
print(xtable(tabISAF))
print("Taliban")
print(xtable(tabTaliban))
sink()

harmApproachCoefs <- read.csv("results/estimatesHarmApproach.ci..individual.fit.csv")

pdf("figs/Fig6-CoefsHarm.pdf", height = 5/3, width = 6.5)
par(mfcol = c(1,3), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.25, 1.75, 3.25, 3.75)


### ISAF Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor ISAF", outer = TRUE, side = 2, cex = 0.65, at = 0.5, line = 1, padj = .1)

coefsPlot <- rbind(indivCoefs[3,], areaCoefs[3,],
                   indivCoefs[1,], areaCoefs[1,]
                   )


points(x.coords, coefsPlot[,"est"], pch = rep(c(19,1), 4), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

text(x.coords[1] , 0.6, "Self or\nfamily", cex = 1)
text(x.coords[2] + .25, -.9, "Manteqa\n(area)", cex = 1)

mtext("Victimization", outer = TRUE, cex = 0.65, at = .482, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = .43, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = .54, line = -0.5)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor Taliban", outer = FALSE, side = 2, cex = 0.65, at = 0, line = 2.8, padj = 0)

coefsPlot <- rbind(indivCoefs[19,], areaCoefs[19,],
                   indivCoefs[17,], areaCoefs[17,]
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,1), 4), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Victimization", outer = TRUE, cex = 0.65, at = 0.16, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.1, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.21, line = -0.5)

### Difference Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Difference in Change\nin Support (Taliban - ISAF)", outer = FALSE, side = 2, padj = -.3, cex = 0.65, at = 0.0, line = 2.5)

coefsPlot <- rbind(indivCoefs[35,], areaCoefs[35,],
                   indivCoefs[33,], areaCoefs[33,]
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,1), 4), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Victimization", outer = TRUE, cex = 0.65, at = 0.82, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.77, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.87, line = -0.5)

dev.off()

###
###
## second plot - approach

indivCoefs <- read.csv("results/estimates.ci..individual.fit.csv")
areaCoefs <- read.csv("results/estimates.ci..area.fit.csv")

harmApproachCoefs <- read.csv("results/estimatesHarmApproach.ci..individual.fit.csv")

pdf("figs/Fig7-CoefsApproach.pdf", height = 5/3, width = 6.5)
par(mfcol = c(1,3), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.33, 1.66, 1.99, 3, 3.33, 3.66 )

### ISAF Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor ISAF", outer = TRUE, side = 2, cex = 0.65, at = 0.5, line = 1, padj = .1)

coefsPlot <- rbind(indivCoefs[3,], ## harm and no approach
                   indivCoefs[7,], ## diff approach - no approach
                    harmApproachCoefs[4,], ## harm and approach

                   indivCoefs[1,], ## harm and no approach
                   indivCoefs[5,], ## diff approach - no approach
                   harmApproachCoefs[1,] ## harm and approach
                   
                   )


points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])


mtext("Experience with Violence", outer = TRUE, cex = 0.65, at = .482, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = .43, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = .54, line = -0.5)

### ISAF Support

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor Taliban", outer = FALSE, side = 2, cex = 0.65, at = 0, line = 2.8, padj = 0)

coefsPlot <- rbind(indivCoefs[19,], ## harm and no approach
                   indivCoefs[23,], ## diff approach - no approach
                   harmApproachCoefs[5,], ## harm and approach

                   indivCoefs[17,], ## harm and no approach
                   indivCoefs[21,], ## diff approach - no approach
                   harmApproachCoefs[2,] ## harm and approach
                   
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Experience with Violence", outer = TRUE, cex = 0.65, at = 0.15, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.1, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.21, line = -0.5)

### Difference Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Difference in Change in\nSupport (Taliban - ISAF)", outer = FALSE, side = 2, padj = -.3, cex = 0.65, at = 0.0, line = 2.5)

coefsPlot <- rbind(indivCoefs[35,], ## harm and no approach
                   indivCoefs[39,], ## diff approach - no approach
                   harmApproachCoefs[6,], ## harm and approach

                   indivCoefs[33,], ## harm and no approach
                   indivCoefs[37,], ## diff approach - no approach
                   harmApproachCoefs[3,] ## harm and approach
                   
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Experience with Violence", outer = TRUE, cex = 0.65, at = 0.82, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.77, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.87, line = -0.5)

par(xpd=NA)
legend(-4.5, -1.4, c("Victimization with No Approach", "Effect of Being Approached", "Victimization and Approach"), 
       pch = c(19,15,17), horiz=T, xjust=0.5, box.lty = 0, pt.cex = 1.3)

dev.off()





###
###
## third plot - pro Taliban

###
## pro taliban plot, 1 x 3
###

indivCoefs <- read.csv("results/estimates.ci..individual.fit.csv")
indivCoefsProTal <- read.csv("results/estimatesInteraction.ci.ProTal.individual.fit.csv")

harmApproachCoefs <- read.csv("results/estimatesHarmApproach.ci.ProTal.individual.fit.csv")

pdf("figs/Fig11-CoefsProTaliban.pdf", height = 5/3, width = 6.5)
par(mfcol = c(1,3), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

##x.coords <- c(1.33, 1.66, 1.99, 3, 3.33, 3.66 )
x.coords <- c(1.25, 1.75, 3.25, 3.75 )

### ISAF Support


### ISAF Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor ISAF", outer = TRUE, side = 2, cex = 0.65, at = 0.5, line = 1, padj = .1)

coefsPlot <- rbind( indivCoefsProTal[91,], ## harm Pro Taliban -- ISAF
                   indivCoefsProTal[3,], ## harm Anti Taliban -- ISAF
                   ##indivCoefsProTal[17,] ## diff pro - anti

                   indivCoefsProTal[88,], ## harm Pro Taliban -- Taliban
                   indivCoefsProTal[1,] ## harm Anti Taliban -- Taliban
                   ##indivCoefsProTal[16,], ## diff pro - anti
                   
                   
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])


mtext("Victimization", outer = TRUE, cex = 0.65, at = .482, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = .43, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = .54, line = -0.5)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Change in Support\nfor Taliban", outer = FALSE, side = 2, cex = 0.65, at = 0, line = 2.8, padj = 0)

coefsPlot <- rbind(
                   indivCoefsProTal[92,], ## harm Pro Taliban -- ISAF
                  indivCoefsProTal[20,], ## harm Anti Taliban -- ISAF
                   ## indivCoefsProTal[34,] ## diff Pro - Anti
                   
                   indivCoefsProTal[89,], ## harm Pro Taliban -- Taliban
                   indivCoefsProTal[18,] ## harm Anti Taliban -- Taliban
                  ## indivCoefsProTal[33,], ## diff Pro - Anti
                   
                   
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Victimization", outer = TRUE, cex = 0.65, at = 0.15, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.1, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.21, line = -0.5)

### Difference Support
plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(0.5, 4.5), ylim = c(-1.2, 1.2), axes = FALSE)

axis(2, at = seq(from = -1, to = 1, by = 0.5), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Difference in Change in\nSupport (Taliban - ISAF)", outer = FALSE, side = 2, padj = -.3, cex = 0.65, at = 0.0, line = 2.5)

coefsPlot <- rbind(indivCoefsProTal[93,], ## harm Pro Taliban -- ISAF
                  indivCoefsProTal[37,], ## harm Anti Taliban -- ISAF
                   ## indivCoefsProTal[51,] ## diff pro - anti

                   indivCoefsProTal[90,], ## harm Pro Taliban -- Taliban
                   indivCoefsProTal[35,] ## harm Anti Taliban -- Taliban
                  ## indivCoefsProTal[50,], ## diff pro - anti
                   
                    
                   )

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15), 2), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

mtext("Victimization", outer = TRUE, cex = 0.65, at = 0.82, line = -0.5, padj = -2)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.77, line = -0.5)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.87, line = -0.5)

par(xpd=NA)
legend(-4, -1.4, c("Pro-Taliban Pashtun", "Non-Pro-Taliban Pashtun"), 
       pch = c(19,15), horiz=T, xjust=0.5, bty = "o", box.lty = "blank", pt.cex = 1.3)

dev.off()

